🎓 Домашняя работа #1 🎓¶


📜 Выполнила:¶

Малкова Ксения Эдуардовна

📧 Email:¶

kemalkova@edu.hse.ru

🗓️ Дата:¶

15 сентября 2024


Файл с заданиями

Замечание: у меня начала заедать буква 'Б' на клавиатуре, поэтому я могла иногда ее пропустить и этого не заметить, ровно как и запятые. Я старалась все тщательно перечитать, но могла не все подправить. Так что заранее извиняюсь за такое написание!

Библиотеки¶

In [1]:
%load_ext autoreload
%autoreload 2
In [2]:
import os
import pandas as pd
import numpy as np
import zipfile
import shutil
from tqdm import tqdm

import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style("darkgrid", {"axes.facecolor": ".9"})

import warnings
warnings.simplefilter('ignore')
In [3]:
from task2.polysomy_search import KlinefelterSearcher
from task3.find_bam_assembly import BamAssemblyFinder
from task4.analyze_gene_cov import GeneCoverageAnalyzer
from task5.find_vcf_assembly import VCFAssemblyFinder
from task6.pair_compare_vcf import VCFComparator
from task7.analyze_genetic_variants import GeneticVariantAnalyzer
from task7.plot_3d_scatter import plot_3d_scatter

Задание 1¶

Определить APOE статус и риск болезни Альцгеймера для полногеномного сиквенса с помощью IGV.

  1. BAM (70Гб)

  2. BAI

Решение¶

  1. Заходим в IGV и выбираем референсный геном hg19
  2. Загружаем BAM и BAI файлы (File $\rightarrow$ Load from URL...)
  3. В IGV в поле поиска вводим ген APOE - таким образом мы перемещаемся в нужную область генома
  4. Определяем аллели APOE. Для этого бежим в SNPедию и понимаем, что основные варианты гена APOE связаны с полиморфизмами rs429358 и rs7412
rs429358 rs7412 Name
C T ε1
T T ε2
T C ε3
C C ε4

Координаты на hg19 ищем в Varsome:

  • rs429358: chr19:45,411,941
  • rs7412: chr19:45,412,079
  1. Возвращаемся в IGV и сравниваем свои данные в этих позициях с референсным геномом:
  • rs429358: на референсном геноме T. У нас гомозиготный вариант T/T
  • rs7412: на референсном геноме C. У нас гомозиготный вариант C/C
  1. Теперь снова идем в SNPедию и оцениваем, какой риск Альцгеймера соответствует нашему генотипу (точнее генотипу нашего образца, а не нашему). Magnitude показывает степень риска, где 6 указывает на очень высокий риск, а 1 - на низкий риск развития АЛьцгеймера.
Common name Genoset Magnitude rs429358 rs7412 Comment
Apo-ε1/ε1 gs267 6 (C;C) (T;T) the rare missing allele
Apo-ε1/ε2 gs271 2.5 (C;T) (T;T)
Apo-ε1/ε3 gs270 2.6 (C;T) (C;T) ambiguous ε2/ε4 or ε1/ε3
Apo-ε2/ε4 gs270 2.6 (C;T) (C;T) ambiguous ε2/ε4 or ε1/ε3
Apo-ε1/ε4 gs272 2.5 (C;C) (C;T)
Apo-ε2/ε2 gs268 4 (T;T) (T;T) good; lowest risk
Apo-ε2/ε3 gs269 2 (T;T) (C;T)
Apo-ε3/ε3 gs246 2 (T;T) (C;C) the most common
Apo-ε3/ε4 gs141 3 (C;T) (C;C)
Apo-ε4/ε4 gs216 6 (C;C) (C;C) ~11x increased Alzheimer's risk

Вывод

Для нашего образца характерен наиболее распространенный статус APOE: Apo-ε3/ε3, который указывает на низкий риск развития Альцгеймера (2 по 6-бальной шкале)

Задание 2¶

  • Найти человека с синдромом Клайнфельтера.
  • Использовать долю прочтений для X, Y хромосомы. Визуализировать образцы на графике.

Data

Решение¶

Идея¶

У мужчин с синдромом Клайнфельтера имеется как минимум одна лишняя X-хромосома (кариотип 47,XXY вместо нормального 46,XY). То есть у мужчин с данным синдромом 2 (или более) X-хромосомы.

У здорового мужчины количество прочтений, приходящихся на X-хромосому, должно быть примерно такое же, как на Y-хромосому. Это означает, что покрытие (coverage) обеих хромосом будет схожим, так как у мужчин по одной X и Y хромосоме. Однако при синдроме Клайнфельтера ожидается, что количество ридов на X-хромосоме будет примерно в 2 раза больше, чем на Y-хромосоме, так как у таких людей две X-хромосомы.

Ход решения¶

  1. Извлечение количества ридов для каждой хромосомы. Для кадого из .bam файлов извлекаем количество ридов, пришедшихся на каждуюиз хромосом. Для этого будем использовать samtools. При этом,

    • $C_X$ - количество прочтений для X-хромосомы
    • $C_Y$ - количество прочтений для Y-хромосомы
  2. Нормализация. Нужно выполнить нормализацию на длину хромосомы. Когда мы говорим о количестве прочтений на участок, то есть сильная зависимость от длины этого участка. Чем он будет длиннее, тем больше ридов на него (скорее всего) придется. Поэтому корректнее будет сравнивать количество прочтений, нормализованное на длину каждой хромосомы. $$N_X = \frac{C_X}{L_X} \quad N_Y = \frac{C_Y}{L_Y},$$ где $L_X$ и $L_Y$ - длина $X$ и $Y$ хромосом соответственно

Ожидаем, что $\text{Ratio} =\frac{N_X}{N_Y}$ будет ранжироваться следующим образом: $$\text{Ratio}[\text{female}] < \text{Ratio}[\text{health\_male}] < \text{Ratio}[\text{disease\_male}]$$

  1. Визуализация. Построим столбчатую диаграмму для каждого образца, которая будет отражать $N_X$ и $N_Y$. Так как у различных образцов различная глубина секвенирования, то для более красивой и корректной визуализации будем нормализовать значения $N_X$ и $N_Y$ на $N_{total}$ (нормализованное суммарное число ридов)
In [4]:
# Указываем путь к папке с .bam и .bam.bai файлами
DATA_FOLDER_TASK2 = "task2/data/downsampled"
In [5]:
searcher = KlinefelterSearcher(data_folder=DATA_FOLDER_TASK2)
data_karyo = searcher.get_normalized_reads(save_data=True)

data_karyo.head(2)
Process .bam files: 100%|██████████| 26/26 [00:00<00:00, 141.78it/s]
Out[5]:
N_X N_Y N_total N_X_norm N_Y_norm Ratio
HG00105 0.000012 6.686477e-06 0.000021 0.556014 0.319458 1.740491
HG00099 0.000024 2.357952e-07 0.000024 1.015929 0.009948 102.124853
In [6]:
searcher.plot_karyotypes(data_karyo, save_fig=True)
No description has been provided for this image

Выводы¶

  • Образцы, $\text{Ratio}$ которых превышает десятки - принадлежат женскому полу (97, 99, 100, 102, 106, 1860). Мы видим какое-то число ридов, которые выровнялись на Y-хромосому, но это случайное событие. Такие риды пришли из других участков генома.
  • Образцы, $\text{Ratio}$ которых меньше двойки - здоровые мужчины (101, 105, 107, 108, 109). Значения $\text{Ratio}$ у них не равняются ровно единице, потому что при нормировке мы не могли учесть все факторы, которые определяют число ридов, которое упадет на тот или иной участок (хромосому). Но с нашим априорным знанием о том, что большинство мужчин здоровы, и результатом в несколько похожих значений $\text{Ratio}$, мы можем заключить, что у здоровых мужчин $\text{Ratio}$ может быть как минимум в диапазоне [1.35, 1.75]
  • Образец 1337 имеет $\text{Ratio} = 3.48$, что отличает его ото всех остальных образцов. Данное значение сильно выше, чем у мужчин (как мы решили, здоровых), но и ниже, чем у женщин. Так что у пациента 1337 кариотип 47,XXY и он обладатель синдрома Кляйнфейтера

Задание 3¶

Определить сборку BAM (на какой из референсных геномов производилось выравнивание)

Data

Решение¶

В принципе в этом задании расписывать особо нечего - просто с помощью samtools достаем метаинформацию о .bam файле (header), где можно найти информацию о сборке файл. Дальше нужно "запарсить" некоторые теги, которые могут содержать информацию о сборке (что указано в документации samtools):

  • @SQ - Reference sequence dictionary, в поле AS=Genome assembly identifier. Сразу скажу, что в этом файле я его не нашла, поэтому:
  • @PG - Program. В общем, можно попытаться извернуться и просмотреть все команды, которыми получали .bam-файл. Там должно быть название референсного генома, в котором может быть отсылка к сборке.

У нас данные в .zip-формате. Чтобы не плодить файлы и не занимать память, извлечем файлы во временную директорию и после удалим, оставив только изначально сжатый формат.

In [7]:
DATA_FOLDER_TASK3 = "task3/data"

task3_bampath = os.path.join(DATA_FOLDER_TASK3, "task3.bam.assembly.zip")
extracted_dir = os.path.join(DATA_FOLDER_TASK3, "tmp")
In [8]:
# 1. Распаковываем архив во временную папку
with zipfile.ZipFile(task3_bampath, 'r') as zip_ref:
    zip_ref.extractall(extracted_dir)

# 2. Достаем из временной папки название .bam файла
for file_name in os.listdir(extracted_dir):
    if file_name.endswith(".bam"):
        bam_file = os.path.join(extracted_dir, file_name)
In [9]:
# Посмотрим на красивое название
bam_file
Out[9]:
'task3/data/tmp/vi0006.genome.markdup.downsampled.bam'
In [10]:
# 3. Наконец-то найдем информацию о сборке
bam_assembly_finder = BamAssemblyFinder(bam_file_path=bam_file)

bam_assembly_finder.find_assembly_info()
Сборка (AS в теге @SQ) не найдена
Найдена информация о программе выравнивания (@PG):
@PG	ID:bwa	PN:bwa	VN:0.7.17-r1188	CL:bwa mem -t 30 -K 10000000 -k 30 -T 30 -R @RG\tSM:vi0006\tID:genome\tPL:illumina /home/xsukhanova/ref/hg38/hg38.fna vi0006.genome.cut_1.fastq.gz vi0006.genome.cut_2.fastq.gz
@PG	ID:samtools	PN:samtools	PP:bwa	VN:1.19.2	CL:samtools view -@ 30 -Sb -
@PG	ID:samtools.1	PN:samtools	PP:samtools	VN:1.19.2	CL:samtools sort -@ 30 -o vi0006.genome.sort.bam -
@PG	ID:MarkDuplicates	VN:Version:4.2.0.0	CL:MarkDuplicates --INPUT vi0006.genome.sort.bam --OUTPUT vi0006.genome.hg38.markdup.bam --METRICS_FILE vi0006.genome.markdup.metrics.txt --MAX_SEQUENCES_FOR_DISK_READ_ENDS_MAP 50000 --MAX_FILE_HANDLES_FOR_READ_ENDS_MAP 8000 --SORTING_COLLECTION_SIZE_RATIO 0.25 --TAG_DUPLICATE_SET_MEMBERS false --REMOVE_SEQUENCING_DUPLICATES false --TAGGING_POLICY DontTag --CLEAR_DT true --DUPLEX_UMI false --ADD_PG_TAG_TO_READS true --REMOVE_DUPLICATES false --ASSUME_SORTED false --DUPLICATE_SCORING_STRATEGY SUM_OF_BASE_QUALITIES --PROGRAM_RECORD_ID MarkDuplicates --PROGRAM_GROUP_NAME MarkDuplicates --READ_NAME_REGEX <optimized capture of last three ':' separated fields as numeric values> --OPTICAL_DUPLICATE_PIXEL_DISTANCE 100 --MAX_OPTICAL_DUPLICATE_SET_SIZE 300000 --VERBOSITY INFO --QUIET false --VALIDATION_STRINGENCY STRICT --COMPRESSION_LEVEL 2 --MAX_RECORDS_IN_RAM 500000 --CREATE_INDEX false --CREATE_MD5_FILE false --GA4GH_CLIENT_SECRETS client_secrets.json --help false --version false --showHidden false --USE_JDK_DEFLATER false --USE_JDK_INFLATER false	PN:MarkDuplicates
@PG	ID:samtools.2	PN:samtools	PP:samtools.1	VN:1.19.2	CL:samtools view -@ 45 -s 0.028 -bo vi0006.genome.hg38.markdup.downsampled.x1.bam vi0006.genome.hg38.markdup.bam
@PG	ID:samtools.3	PN:samtools	PP:samtools.2	VN:1.20	CL:samtools view -h -b -o vi0006.genome.markdup.downsampled.bam vi0006.genome.markdup.downsampled.x1.bam chr20
@PG	ID:samtools.4	PN:samtools	PP:samtools.3	VN:1.20	CL:samtools view -H task3/data/tmp/vi0006.genome.markdup.downsampled.bam
In [11]:
# 4. Удаляем временную папку
shutil.rmtree(extracted_dir)

Вывод¶

Ура! Мы нашли полную команду запуска выравнивания, где было явно указано, что геном сборки hg38!

Задание 4¶

Из BAM файла получить BED файл с регионами, покрытие которых 10 и более. С помощью bedtools определить, какая доля гена APOE имеет покрытие x10+.

Data

Решение¶

  1. Получение информации о покрытии и фильтрация - с помощью bedtools genomecov получаем .bed файл с покрытием регионов, после чего отбираем только те, где покрытие >x10.
  2. Определение координат APOE
    • Определение сборки - с помощью samtools в .bam файле находим информацию о сборке генома
    • Поиск координат APOE в соответствующей сборке генома - можно было просто зайти на NCBI, или его запарсить. Но! Есть прекрасная библиотека в питоне mygene, которая делает все сама. Вообще MyGene — это веб-тул, который построен поверх множества баз данных (NCBI, Ensembl, UniProt, ClinVar etc). Через REST API он отправляет запросы в эти БД по гену (или его символу, или NCBI/RefSeq идентификатору) и подтягивает нужную нам информацию (координаты для разных сборок, названия гена, варианты и связь с заболеваниями (в случае ClinVar и подобных БД), Gene Ontology и прочее). В общем, крутая штука! Горжусь братьями по делу!
  3. Пересечение APOE и найденных регионов с покрытием >x10
    • Создание .bed файла с координатами APOE
    • Поиск пересечений - с помощью bedtools intersect вытаскиваем участки из нашего .bed файла, соответствующие APOE
In [12]:
DATA_FOLDER_TASK4 = "task4/data"

task4_bampath = os.path.join(DATA_FOLDER_TASK4, "task4.bam")

Воспользуемся решением прошлого задания для определения сборки .bam-файла

In [13]:
bam_assembly_finder = BamAssemblyFinder(bam_file_path=task4_bampath)
_ = bam_assembly_finder.find_assembly()
Найдена сборка (по тегу @SQ в поле AS): NCBI37

Странная аннотация, но окей, похоже на GRCh37. Так что кажется наша сборка - hg19

P.S.: здесь тоже удивились такому странному id, потому что он не указывает, что это именно человека :)

Теперь наконец-то найдем долю гена APOE с покрытием >x10. Для этого воспользуемся GeneCoverageAnalyzer (подробнее в файле task4/analyze_gene_cov.py)

In [14]:
gene_cov_analyzer = GeneCoverageAnalyzer(bam_file_path=task4_bampath, 
                                         gene="APOE",
                                         coverage_threshold=10,
                                         genome_assembly="hg19",
                                         species="human")

gene_cov_info = gene_cov_analyzer.compute(return_info=True)
Файл с координатами гена APOE уже существует по пути task4/data/APOE_human.bed
Файл с отфильтрованными регионами по покрытию уже существует по пути task4/data/task4_x10.bed
Файл пересечения уже существует по пути task4/data/task4_APOE_human.bed
Длина гена APOE: 3639 п.о.
Длина пересечения гена APOE и .bam-файла: 647 п.о.
Доля гена APOE с покрытием больше 10x: 17.78%
In [15]:
gene_cov_info
Out[15]:
{'intersected_length': 647,
 'gene_length': 3639,
 'coverage_fraction': 0.17779609782907393}

Я сразу оговорюсь - здесь доля покрытия может чуть варьироваться, потому что:

  1. Покрытие больше x10 - включая или нет, не очень ясно, я брала больше или равно 10x (не включая 10x было на уровне 14.5 %)
  2. На Ensembl и NCBI) чуть разные точки старта APOE в сборке hg19. В моей версии Mygene() подтянул вариант из Ensembl
In [16]:
!cat {gene_cov_analyzer.gene_coord_file_path}
19	45409011	45412650

Вывод¶

Какое-то достаточно низкое покрытие.

Задание 5¶

Определить сборку (референсный геном) VCF файла.

Data

Решение¶

В принципе схема аналогичная, как в задании №3. Открываем header .vcf, только уже с помощью bcftools, и ищем там информацию о сборке. Можно двумя путями:

  • Ищем assembly напрямую (я видела это в информации о контигах)
  • Лио ищем референс - навание файла, в котором может быть отражена сборка
In [17]:
DATA_FOLDER_TASK5 = "task5/data"

task5_vcfpath = os.path.join(DATA_FOLDER_TASK5, "task5.input.vcf.assembly.zip")
extracted_dir = os.path.join(DATA_FOLDER_TASK5, "tmp")

Эх и снова .zip. Делаем ту же схему, что и выше: распаковка $\rightarrow$ временная папка $\rightarrow$ делаем дела $\rightarrow$ удаляем временную папку

In [18]:
# Распаковка файла во временную папку
with zipfile.ZipFile(task5_vcfpath, 'r') as zip_ref:
    zip_ref.extractall(extracted_dir)
    print(f'Files: {os.listdir(extracted_dir)}')
Files: ['ki3640.imputed.vcf.gz']
In [19]:
# Указываем путь к файлу, находящимся во временной папке
task5_vcfpath = os.path.join(extracted_dir, "ki3640.imputed.vcf.gz")
In [20]:
# 3. Наконец-то найдем информацию о сборке
vcf_assembly_finder = VCFAssemblyFinder(vcf_file_path=task5_vcfpath)

vcf_assembly_finder.find_assembly_info()
Найдена сборка (по тегу assembly): hg19

И снова hg19!

Можем попробовать и по референсному файлу, но вообще это вариант, когда вариант с поиском assembly не сработает:

In [21]:
vcf_assembly_finder.find_reference()
Найдено референс (по тегу ##reference): file:///mnt/db/hg19/hg19.fa
##source="beagle.18May20.d20.jar"
##FILTER=<ID=LowQual,Description="Low quality"
Out[21]:
'file:///mnt/db/hg19/hg19.fa\n##source="beagle.18May20.d20.jar"\n##FILTER=<ID=LowQual,Description="Low quality"'

Не забываем удалять временную папку:

In [22]:
shutil.rmtree(extracted_dir)

Задание 6¶

Попарно сравнить генотипы для 3 образцов с помощью bcftools. Сделать выводы о родстве.

Data

Решение¶

Для попарного сравнения генотипов будем использовать bcftools gtcheck. И здесь встает важный вопрос - запускаем с флагом -E 0 или без него?

  1. По умолчанию работаем с флагом -E 40:
    • Существует ненулевая ($10^{-4}$) вероятность того, что несовпадение между генотипами вызвано ошибкой генотипирования.
    • В таком случае Discordance может быть и дробным числом
  2. С флагом -E 0:
    • Каждое несоответствие между генотипами считается жестким несовпадением.
    • В таком случае Discordance может быть инетрпретирован как количество несовпадающих генотипов (ссылка), что число всегда целое.

Этому вопросу посвящено очень много обсуждений на гитхабе (#1538, #1556, #2194). И люди радовались, когда получали точное число несовпадающих генотипов (хотя отношение между Discordance для разных пар сохраняется). В общем, используем флаг -E 0*!

*Раньше -E был -e. В новой версии bcftools -e стал означать --exclude, так что нужно использовать тот или иной флаг в зависимости от своей версии bcftools.

**Кстати есть прикол, что bcftools gtcheck работает только с биаллельными сайтами, а мультиаллельные просто игнорирует (что обсуждалось в #1152 и сообщено в их коде). Можно предварительно "разбить" все мультиаллельные сайты на несколько строчек (bcftools norm -m-), как бы сделая их биаллельными. Я это прогоняла, и в наших образцах в принципе нет мультиаллельных, так что в нашем случае эта процедура избыточна :)

In [23]:
DATA_FOLDER_TASK6 = "task6/data"
task6_vcfpath = os.path.join(DATA_FOLDER_TASK6, "task6.vcf.gz")
In [24]:
comparator = VCFComparator(vcf_file_path=task6_vcfpath)
variant_comparison_df =  comparator.compare(save_result=True)
Результат сравнения образцов из .vcf записан в task6/data/pairwise_comparisons.csv
In [25]:
variant_comparison_df
Out[25]:
Sample1 Sample2 Discordance Average_neg_log_P_HWE Number_of_sites_compared Number_of_matching_genotypes
0 HG00732 HG00733 41781.0 0.030539 1066557 1024776
1 NA20320 HG00733 73481.0 0.014936 1066557 993076
2 NA20320 HG00732 74732.0 0.014677 1066557 991825

Что эти колонки вообще обозначают?¶

  1. Sample1

Первый образец (его идентификатор) из сравниваемой пары

  1. Sample2

Второй образец (его идентификатор) из той же пары, с которым сравнивается первый образец

  1. Discordance

Метрика, отражающая, насколько сильно отличаются генотипы в паре образцов. В нашем случае, а именно при использовании флага -E 0, Discordance - это количество несовпадающих генотипов.

  1. Average_neg_log_P_HWE

Вероятность отклонения (а точнее минус логарифм от нее) генотипов от равновесия Харди-Вайнберга (HWE). Согласно HWE, в идеальной популяции, где спаривание происходит случайным образом и отсутствуют такие факторы как мутация, миграция, отбор и генетический дрейф, частоты генотипов остаются стабильными от поколения к поколению.

Хотя эти предположения редко бывают верны в естественном мире, они позволяют нам вычислить ожидаемую частоту аллеля. Значительные различия между наблюдаемыми и ожидаемыми частотами указывают на то, что «что-то» (т. е. одно или несколько из вышеперечисленного) происходит .

В контексте результата bcftools, для каждой позиции (локуса) в геноме у нас есть частота каждого аллеля и частота генотипов в этой позиции. На основе этих частот можно рассчитать ожидаемые частоты генотипов при условии HWE. Затем наблюдаемые частоты сравниваются с ожидаемыми, и на основе этого сравнения для каждой позиции оценивается P(HWE) - вероятность того, что частоты генотипов, наблюдаемые для этого сайта, соответствуют равновесию Харди-Вайнберга.

В bcftools, когда выполняется проверка на соответствие равновесию Харди-Вайнберга, используется $\chi^2$-тест:

$$\chi^2 = \sum{\frac{(O_i - E_i)^2}{E_i}}$$

где $O_i$ - наблюдаемые частоты, а $E_i$ - ожидаемые частоты

  1. Number_of_sites_compared

Общее количество генетических позиций (сайтов), по которым было выполнено сравнение между двумя образцами

  1. Number_of_matching_genotypes

Количество позиций, по которым генотипы образцов совпадают - сколько генетических позиций идентичны между двумя сравниваемыми образцами. В случае нашего флага -E 0: $$\text{[\text{Number of sites compared}]} = \text{[Number of matching genotypes]} + \text{[Discordance]}$$

На этом моменте я искренне перестала понимать, почему люди с обсуждений на гите так счастливы использовать флаг -E 0, если количество несовпадающих генотипов можно и самому рассчитать... а с тем же дефолтным -E 40 можно хотя бы регулировать вероятность ошибки генотипирования

In [26]:
variant_comparison_df['Discordance'] + \
    variant_comparison_df['Number_of_matching_genotypes'] \
        == variant_comparison_df['Number_of_sites_compared']
Out[26]:
0    True
1    True
2    True
dtype: bool

Визуализация результата¶

1. Относительное число несоответствий¶

Нарисуем тепловую карту (heatmap) по Discordance. Более "яркие" значения соответствуют большим значениям Discordance - то есть такая пара образцов имеет большое количество несовпадающих генотипов, и наоборот. Самый маленький Discordance у пары образцов сам-с-собой (sample vs sample) и будет равным 0

In [27]:
# Получаем уникальные пары оразцов
samples = pd.unique(variant_comparison_df[['Sample1', 'Sample2']].values.ravel('K'))

# Создаем пустую симметричную матрицу с нулями
matrix = pd.DataFrame(np.zeros((len(samples), len(samples))), index=samples, columns=samples)

# Заполняем матрицу значениями Discordance
for _, row in variant_comparison_df.iterrows():
    matrix.at[row['Sample1'], row['Sample2']] = row['Discordance']
    matrix.at[row['Sample2'], row['Sample1']] = row['Discordance']

# И еще нормализуем на общее число сравниваемых сайтов - оно одинаково между образцами
# Идея такая, что хочется посмотреть на % несовпадающих сайтов, а не их абсолютное число
assert variant_comparison_df['Number_of_sites_compared'].nunique() == 1, \
    "'Number_of_sites_compared' различается между образцами"

n_sites = variant_comparison_df['Number_of_sites_compared'].iloc[0].item()
matrix = 100 * matrix / n_sites
In [28]:
ax = sns.heatmap(matrix, 
                 cmap = 'Purples', linewidths=1, linecolor="grey",
                 annot=True, fmt=".2f")
ax.set_xticklabels(matrix.columns, rotation=45)
ax.set_yticklabels(matrix.index, rotation=0)
ax.set_title("Относительное число несовпадений (%)")

plt.show()
No description has been provided for this image

Самое большое число несоответствий у пары NA20320-HG00732, и примерно такое же число несоответствий - у пары NA20320-HG00733. Конечно наша выборка только из трех образцов (и $C_3^2=3$ пар образцов), но что-то мне подсказывает, что NA20320 сильно отбивается от HG00732 и HG00733. Точнее не что-то, а 1) названия образцов, 2) значения Discordance

2. Соблюдение HWE¶

Не будем отклоняться от прошлого графика - также нарисуем тепловую карту средних значений вероятности отклонения от HWE

In [29]:
# Создаем пустую симметричную матрицу для Average_neg_log_P_HWE
hwe_matrix = pd.DataFrame(np.zeros((len(samples), len(samples))), index=samples, columns=samples)

# Заполняем матрицу значениями Average_neg_log_P_HWE
for _, row in variant_comparison_df.iterrows():
    hwe_matrix.at[row['Sample1'], row['Sample2']] = row['Average_neg_log_P_HWE']
    hwe_matrix.at[row['Sample2'], row['Sample1']] = row['Average_neg_log_P_HWE']
In [30]:
# Визуализируем тепловую карту
ax = sns.heatmap(hwe_matrix, 
                 cmap='Blues', linewidths=1, linecolor="grey",
                 annot=True, fmt=".2f")
ax.set_xticklabels(hwe_matrix.columns, rotation=45)
ax.set_yticklabels(hwe_matrix.index, rotation=0)
ax.set_title(r"AVG $\left[ -\log{P(\text{HWE})} \right]$", fontsize=14)

plt.show()
No description has been provided for this image

Чтобы интерпретировать картинку, сделаем небольшое вступление (на всякий случай, я сама порой путаюсь):

  • $P(\text{HWE})$ - вероятность отклонения от HWE для пары образцов.

    • Бóльшая вероятность указывает на соблюдение HWE у двух образцов по этой позиции
    • Мéньшая вероятность указывает на отклонение от HWE у двух образцов по этой позиции
  • $-\log{P(\text{HWE})}$ - то же самое, что и выше, только с мат. наворочками над исходной вероятностью:

    • Низкие значения соответствуют более высокому p-value, что указывает на то, что нет значимых отклонений от HWE у двух образцов по этой позиции.
    • Высокие значения соответствуют более маленькому p-value, что указывает на то, что с какой-то вероятностью можно говорить о значимых отклонениях от HWE в этой паре образцов по данной позиции.

На нашей тепловой карте мы видим только низкие значения -$\log{P(\text{HWE})}$, что указывает на отсутствие значительных отклонений от равновесия Харди-Вайнберга среди всех пар образцов. А если у нас нет значительных нарушений HWE, то это может указывать на отсутствие значительных генетических различий или мутаций между образцами в сравниваемых позициях.

Выводы¶

Мы посмотрели на красивые картинки и сделали два отдельных вывода:

  • Для пар с образцом NA20320 относительное число несовпадений составляет примерно 7%, тогда как для последней пары — около 4%. Эти значения могут показаться значительными на фоне числа сравниваемых позиций, которое составляет порядка миллиона
  • Все пары образцов показывают низкие значения $-\log{P(\text{HWE})}$, что указывает на незначительные отклонения от равновесия HW в рассматриваемых позициях

Таким образом, несмотря на высокое абсолютное и относительное количество несовпадающих позиций между образцами, низкие значения $-\log{P(\text{HWE})}$ свидетельствуют о том, что эти несовпадения не являются статистически значимыми. То есть различия в генотипах не указывают на значимые генетические изменения или мутации, а скорее связаны с естественным генетическим разнообразием.

В общем, данные свидетельствуют о наличии генетического разнообразия между образцами, но не о значительных мутациях, что может указывать на близкое или умеренное родство между ними.

Задание 7¶

Рассчитать частоты генетических вариантов для трех популяций (AFR, EUR, SAS) и построить scatterplot попарно между популяциями (ось Х – частота в первой популяции, ось Y – частота во второй популяции, точка – генетический вариант с соответствующими частотами)

Data

Решение¶

In [31]:
DATA_FOLDER_TASK7 = "task7/data"

task7_path = os.path.join(DATA_FOLDER_TASK7, "task7.zip")
extracted_dir = os.path.join(DATA_FOLDER_TASK7, "tmp")

Мы снова имеем дело с .zip-архивом, поэтому следуем той же логике, что и ранее: распакуем его во временную папку, сделаем дело, и после временную папку подчистим. Я написала кастом для сохранения таблицы частотности аллелей, поэтому результат можем сохранить. Но промежуточное точно не нужно

In [32]:
with zipfile.ZipFile(task7_path, 'r') as zip_ref:
    zip_ref.extractall(extracted_dir)

    print(f'Files: {os.listdir(extracted_dir)}')
Files: ['data']

О боже, папка в папке.. Ну ладно, сделаем ход конем и вручную переместим файлы в tmp/

In [33]:
!mv {os.path.join(extracted_dir, "data/*")} {extracted_dir}
!rm -r {os.path.join(extracted_dir, "data/")}

Супер! Теперь, наконец, определяем пути к файлам. Только перед этим посмотрим на них:

In [34]:
os.listdir(extracted_dir)
Out[34]:
['ALL.chr15.subset.hg38.annotated.clean.vcf.gz', 'tubes_pop.tsv']
In [35]:
population_data_path = os.path.join(extracted_dir, 'tubes_pop.tsv')
frequency_data_path = os.path.join(extracted_dir, 'ALL.chr15.subset.hg38.annotated.clean.vcf.gz')
af_savepath = os.path.join(DATA_FOLDER_TASK7, 'AF_populations.csv')

И приступим к самому главному - подсчитаем частоты! (под капотом в основном bcftools)

In [36]:
genvar_analyzer = GeneticVariantAnalyzer(population_data_path=population_data_path,
                                         frequency_data_path=frequency_data_path)
Количество образцов в анализируемых популяциях:
pop
AFR    671
EUR    521
SAS    492
Name: count, dtype: int64
In [37]:
AF_frequencies = genvar_analyzer.calculate_af(save_path=af_savepath)
Подсчет генетических частот...: 100%|██████████| 3/3 [05:35<00:00, 111.67s/it]
In [38]:
# посмотрим, как выглядит итоговая табличка
AF_frequencies.head(5)
Out[38]:
CHROM POS REF ALT AF Population
0 15.0 17009562.0 T C 0.000000 EUR
1 15.0 17009590.0 G T 0.445298 EUR
2 15.0 17009613.0 A C 0.000000 EUR
3 15.0 17009614.0 T C 0.000960 EUR
4 15.0 17009622.0 G A 0.000000 EUR

Перед тем, как мы предпримем попытку что-то нарисовать и интерпретировать, посмотрим, что у нас вообще за популяции. И здесь интересно, зачем образцы из SGDP отделили... возможно, чтобы контролировать батч-эффект, но это не точно

In [39]:
genvar_analyzer.population_data[['pop', 'pop_name']].value_counts()
Out[39]:
pop  pop_name                              
AFR  African Ancestry                          663
EUR  European Ancestry                         514
SAS  South Asian Ancestry                      486
AFR  African Ancestry,Africa (SGDP)              8
EUR  European Ancestry,West Eurasia (SGDP)       7
SAS  South Asia (SGDP),South Asian Ancestry      6
Name: count, dtype: int64

График рассеяния (scatter plot)¶

In [40]:
genvar_analyzer.plot_scatter(df_frequencies=AF_frequencies, populations=genvar_analyzer.populations)
No description has been provided for this image

красота 👍

Что за точки и как на это смотреть?¶

Думаю не будет вредно разобраться, что означает каждое положение точки на графике. Из формулировки задания освежим, что точка - это генетический вариант с соответствующими частотами (частоты варианта в популяции, отраженной на осях OX и OY). Цветом обозначена позиция (чем темнее - тем ближе к началу 15-й хромосомы, чем светлее - тем ближе к концу 15-й хромосомы). Но так как вариантов очень много, последние точки (т.е. расположенные ближе к концу 15-й хромосомы, и, соответственно более желтые) легли на графике на все остальные. Тем не менее по краям мы можем увидеть точки, соответствующие середине или началу хромосомы :)

Итак, к графику (этот текст я писала, когда я неправильно доставала генетические частоты, и график был более красивый, но сейчас тоже +- актуально):

  1. Точка лежит на линии регрессии или близка к ней, т.е. $y=x$. Если точка находится на этой линии, то в обеих популяциях ~одинаковая распространненность данного генетического варианта. То есть вариант:

    • Возможно имеет стабильное распределение вне зависимости от популяции
    • Возможно нейтральная мутация, не подвергающаяся селекционному давлению ни в одной из популяций внутри пары.
  2. Точка сильно выше линии $y=x$. Тогда генетический вариант более распространен во второй популяции (OY) по сравнению с первой (OX). То есть вариант:

    • Возможно более благоприятен в среде второй популяции, где он дает преимущество (в то время как во второй популяции данный вариант или вреден, или нейтрален $\rightarrow$ меньше частота варианта в этой популяции)
    • Скорее всего вариант имеет важное биологическое значение и подвергается разному давлению выбора в этой паре популяций
  3. Точка сильно ниже линии $y=x$. Тогда генетический вариант более распространен в первой популяции (OX) по сравнению со второй (OY). В принципе схема та же, что и в п.2

Выводы по графику рассеяния¶

  • Если бы мы построили график частот аллелей в одной популяции против частот тех же аллелей в этой же популяции, мы бы получили прямую линию ($y=x$). Соответственно, графики, которые больше прижаты к линии $y=x$, указывают на более высокую генетическую схожесть между сравниваемыми популяциями
  • График рассеяния для пары популяций SAS-EUR наиболее близок к линии $y=x$. Так что, вероятно, популяции SAS и EUR являются наиболее генетически схожими из всех (трех) исследуемых пар популяций
  • Графики для пар популяций AFR-EUR и AFR-SAS занимают значительную часть графика, почти заполняя квадрат [0, 0] - [1, 1]. То есть между этими парами популяций существует значительное различие в частотах генетических вариантов. При этом в обеих парах наблюдается как высокая частота аллелей в популяции AFR и низкая в другой (SAS/EUR), так и наоборот. Скорее всего популяция AFR значительно отличается от популяций EUR и SAS, возможно из-за сильных различий в среде обитания

3D график рассеяния¶

Мне просто стало интересно...

Если не запускаете ноутбук (который, возможно, и не запустится из-за зависимостей, которые на другой машине не могу проверить + некоторые файлы я переименовывала), то можете открыть и поиграться с интерактивными картинками, которые я сохранила в формате .html. Они находятся по этой ссылке

Дальше есть два варианта действий:

  1. Скачать любой файл и открыть локально. Но не факт, что откроется с телефона
  2. Добавить приставку https://html-preview.github.io/?url= к ссылке на .html` Например: https://html-preview.github.io/?url=https://github.com/ksumarshmallow/Practical-BF/blob/main/HW1/task7/imgs/AF_3D_17009562_25506717.html
In [41]:
# Собираем частотные таблицы по каждой популяции
dfs = genvar_analyzer._split_and_filter_dfs(df_frequencies=AF_frequencies, populations=genvar_analyzer.populations)

# Названия популяций
pops = [k for k in dfs.keys()]

# AF для каждой популяции
AFs = np.array([df.AF.values for df in dfs.values()])

# Позиции
positions = dfs[pops[0]].POS.values.astype(int)

AFs.shape
Out[41]:
(3, 2294447)

Все позиции мы не сможем нарисовать (все сломается и будет некрасиво), поэтому нарисуем несколько графиков!

In [42]:
# Создаем папку для сохранения .html, если ее еще нет
AF_3d_folder = os.path.join('task7', 'imgs')

if not os.path.exists(AF_3d_folder):
    os.makedirs(AF_3d_folder)

Наконец-то рисуем интерактив! Будем идти с шагом 100.000 по генетическим вариантам

In [43]:
step = 100_000
kwargs = {"x_title": pops[0], "y_title": pops[1], "z_title": pops[2]}

for i in tqdm(range(0, len(positions), step), desc="Plot 3D AFs...", colour="CYAN"):
    x = AFs[0][i:i+step]
    y = AFs[1][i:i+step]
    z = AFs[2][i:i+step]
    pos = positions[i:i+step]
    plot_name = f'AF_3D_{pos[0]}_{pos[-1]}'

    # рисуем график
    fig = plot_3d_scatter(x=x, y=y, z=z, c=pos, plot_name= plot_name, **kwargs)

    # сохраняем в папочку
    fig.write_html(f"{AF_3d_folder}/{plot_name}.html")
Plot 3D AFs...: 100%|██████████| 23/23 [00:03<00:00,  6.60it/s]

Посмотрим на последнюю картинку, делая сначала ноутук интерактивным:

  • В другой среде возможно нужно будет поставить %matplotlib widget
  • Для просмотра с гитхаба инструкция выше
In [44]:
%matplotlib notebook
fig

Особенности на 3D графике рассеяния¶

  • Лучше на него смотреть, ориентивовав "на себя" главную диагональ
  • Можно регулировать график - наверху для этого есть виджеты (Zoom, Pan=перемещать сам кубик, Orbital rotation, Turntable rotation)
  • Если открыли в интерактиве - можно навести курсор на точку и посмотреть, какие частоты соответствуют этому варианту у всех трех популяций
  • Отклонения от диагонали - варианты, которые отличаются как минимум между двумя популяциями. Они нам наиболее интересны!

Выбросы (outliers)¶

Меня очень сильно понесло (ночное вдохновение) и сейчас я хочу посмотреть на "выбросы" - те генетические варианты, которые отбиваются от остальной кучи точек на графике рассеяния и 3D графике рассеяния.

В целом выбросы интересны потому, что они являются вариантами, где частота аллеля сильно отличается между популяциями. То есть это что-то индивидуальное у одной популяции и может быть следствием особых условий среды в этой популяции. Вообще мне легче будет найти выбросы из 3D граффике рассеяния. Поэтому формулирую следующее:

Цель данного раздела - найти генетические варианты, частота которых различается хотя бы между двумя популяциями (из EUR, AFR, SAS).


В контексте выбросов по генетическим вариантам самый популярный пример - вариант, связанный с геном HBB, который кодирует β-цепь гемоглобина. Мутации в последнем:

  1. Приводят к различным заболеваниям крови, включая серповидноклеточную анемию;
  2. Оказывает сильное влияние на восприимчивость к малярии.

То есть, существует связь между серповидно-клеточной анемей и защитой от малярии. Этот вариант особенно распространен среди жителей Западной и Центральной Африки, некоторых районов Индии и Средиземноморья, где малярия была исторически распространена. Люди с одной копией мутантного аллеля (гетерозиготы) имеют частичную защиту от малярии, что делает этот вариант частым в популяциях, подвергшихся сильному давлению естественного отбора.

Так что это типичный пример того, как различия в окружающей среде могут приводить к резкому расхождению частот аллелей между популяциями. И я хочу посмотреть, что же происходит у нас!


Выбросы заметны на построенном 3D scatter-plot'е, но их нам нужно как-то вычислить. Я пойду наивным образом:

  1. Извлечем частоты для каждого генетического варианта в каждой популяции
  2. Вычислим дисперсию генетических частот среди популяций (да, по трем, но а что делать)
  3. Генетические варианты с самой высокой дисперсией - те самые варианты, частоты которых сильно различаются между тремя популяциями
In [45]:
# Первый шаг - сделан выше

# Второй шаг - считаем стандартное отклонение
deviations = AFs.std(axis=0)
print(f'Максимальная дисперсия: {deviations.max():.4f}')

# Третий шаг:
# 1. Cчитаем порог по отклонению
dev_threshold = deviations.max() - deviations.std() # np.quantile(deviations, 0.995).item()
print(f'Порог по дисперсии: {dev_threshold:.4f}')

# 2. Определяем позиции, где отклонение > порога
deviations_filtered_idx = np.where(deviations > dev_threshold)[0]
positions_filtered = positions[deviations_filtered_idx]
print(f'Количество генетических вариантов после отсечки: {positions_filtered.shape[0]}')

# 3. Находим в них соответствующие частоты
AFs_filtered = AFs[:, deviations_filtered_idx]
Максимальная дисперсия: 0.4175
Порог по дисперсии: 0.3847
Количество генетических вариантов после отсечки: 11

Выведем наши самые крутые позиции!

In [46]:
outliers_variants = pd.DataFrame(AFs_filtered, index = pops, columns = [f'AF_{i+1}' for i in range(len(positions_filtered))]).T
outliers_variants['chr15_pos'] = positions_filtered
outliers_variants
Out[46]:
EUR SAS AFR chr15_pos
AF_1 0.050864 0.019309 0.850969 25275198
AF_2 0.048944 0.019309 0.851714 25275760
AF_3 0.008637 0.001016 0.890462 29135197
AF_4 0.972169 0.938008 0.119225 33966633
AF_5 0.972169 0.938008 0.131893 33967969
AF_6 0.976008 0.940041 0.137854 33973629
AF_7 0.976967 0.940041 0.142325 33977774
AF_8 0.976967 0.958333 0.151267 33979447
AF_9 0.935701 0.966463 0.102832 45095352
AF_10 0.926104 0.960366 0.093890 45105368
AF_11 0.077735 0.115854 0.913562 60563356

Ладно, мой способ нахождения был не самым лучшим, потому чтос самая высокая дисперсия будет там, где хотя бы в паре популяций у одной частота варианта будет >0.9, а у второй <0.1. Но все равно мы нашли, по-моему, что-то интересное! То есть именно те варианты, которые часто встречаются в одной популяции, но редко встречаются в других!. Либо наоборот - редко встречаются в одной популяции, но часто бывают в других. И, кстати, снова от других популяций отличается только AFR

Теперь мы бежим в NCBI Genome Data Viewer и смотрим на наши позиции - возможно, для них уже что-то известно. Для всех нашлись rsid. И табличка по AF среди популяций в любых случаях совпадает с нашей!

Я все обобщила тут:

Название гена Позиция rsID и ссылка Тип (SNV/индел) Публикации ClinVar
SNHG14 chr15:25275198 rs1402074614, rs1442810789, rs2073056425, rs577002978 Индел, SNV Нет публикаций Нет
chr15:25275760 rs2714786 SNV Нет публикаций Нет
chr15:29135197 rs10152250 SNV Нет публикаций Нет
chr15:33966633 rs600637 SNV Нет публикаций Нет
AVEN, CHRM5 2KB Upstream Variant chr15:33967969 rs676281 SNV Нет публикаций Нет
chr15:33973629 rs1175001719, rs476814 Индел, SNV Нет публикаций Нет
chr15:33977774 rs480616 SNV Нет публикаций Нет
chr15:33979447 rs592585 SNV Нет публикаций Нет
DUOX2 chr15:45095352 rs199138 SNV 1. A Population-Based Genomic Study of Inherited Metabolic Diseases Detected Through Newborn Screening
2. Effective selection of informative SNPs and classification on the HapMap genotype data
ClinVar (rs199138)
chr15:45105368 rs269863 SNV Нет публикаций ClinVar (rs269863)
RORA, RORA-AS1 chr15:60563356 rs974828 SNV Нет публикаций Нет

Очищаем пространство¶

In [47]:
# Удаляем временную папку
shutil.rmtree(extracted_dir)

Вывод¶

  1. Мы научились вычислять частоты генетических вариантов для различных популяций
  2. Мы визуализировали частоты на графике рассеяния и 3D графике рассеяния (и поигрались с последним)
  3. Из визуализации частот мы поняли, что популяция AFR отличается по частоте генетических вариантов от популяций EUR и SAS
  4. Мы нашли наиболее варьирующиеся позиции, которые я назвала выбросами - те генетические варианты, которые стоят далеко от главной диагонали на 3D графике рассеяния.
  5. Каждую позицию мы просмотрели в NCBI DATA Viewer, определив, в каком гене находится данный генетический вариант, и найдя соответствующий rsid
  6. Еще раз убедились, что популяция AFR более всего отличается от других двух популяций

Итог¶

Я исчерпала все свое вдохновение и силы программировать, но я в восторге! 🌟